function [Price, n, V] = americanPen( r, sigma, S0, K, T, steps, flag, NodeM, D, dnorm, stepMethod, tol)
% euVarTime Summary of this function goes here
% This function uses variable timestepping to calculate PDE of European
% Option prices. Method used is Crank-Nicolson Timestepping
% with either constant timestepping or variable timestepping
% 
% steps: if constant timestepping is used, steps determines dt
%        if variable timestepping, steps determines starting dTao
% flag = 0: put (Default), flag = 1: call
% NodeM: Node multiplier, default is 1
% stepMethod: 0: constant timesteps; 1: variable timesteps;
% tol: tolerance level, default is 10^-6;
% Returned Value: n: if variable steps, time step tracker

% Argument check
if nargin < 7
    flag = 0; NodeM = 1; D = 1; dnorm = 0.1; stepMethod =0; tol = 10^-6; end
if nargin < 8
    NodeM = 1; D = 1; dnorm = 0.1; stepMethod =0; tol = 10^-6;  end
if nargin < 9
    D = 1; dnorm = 0.1; stepMethod =0; tol = 10^-6; end
if nargin < 10
    dnorm = 0.1; stepMethod =0; tol = 10^-6; end
if nargin < 11
    stepMethod =0; tol = 10^-6; end
if nargin < 12
    tol = 10^-6; end

% Generate S grid
S = [0:0.1*K:0.4*K, 0.45*K:0.05*K:0.8*K, ...
    0.82*K:0.02*K:0.9*K, 0.91*K:0.01*K:1.1*K, ...
    1.12*K:0.02*K:1.2*K, 1.25*K:0.05*K:1.6*K, ...
    1.7*K:0.1*K:2*K, 2.2*K, 2.4*K, 2.8*K, ...
    3.6*K, 5*K, 7.5*K, 10*K];
% Delta Tao
dt = T/steps;
% Number of Asset Prices Grid = Q
Q = length(S);

% Now if we want to increase size of S grid
if NodeM ~= 1
    for i = 2:NodeM
        S_old = S;
        Q = length(S_old);
        S_new = zeros(1,2*Q-1);
        S_new(1:2:end) = S_old;
        S_new(2:2:end-1) = (S_old(1:end-1) + S_old(2:end))/2;
        S = S_new;
        Q = length(S);
    end
end

% Now get the alpha and beta vector
[Alpha, Beta] = alphaBeta(S, r, sigma);

% Now create the sparse matrix M
% This time we use the preallocated indexing sets to create M
% Which is more efficient
i = 2:Q-1;
M = sparse(i,i, (Alpha + Beta + r), Q, Q, 3*Q) + ...
    sparse(i,i-1, -Alpha, Q, Q, 3*Q) + ...
    sparse(i,i+1, -Beta, Q, Q, 3*Q);
M(1,1) = r;

% Option Value at tao = T-T = 0, Default for put option
V = zeros(Q,1);
if flag == 0
    V(Q) = 0;
    V(1:end-1) = max(K - S(1:end-1),0);
else
    V(Q) = S(end);
    V(1:end-1) = max(S(1:end-1) - K, 0);
end
% Now iterate through time steps to solve Vn
Vs = V; % payoff of early exercise V*
if stepMethod == 0 % constant timesteppting
    Mhat = M*dt;
    for i = 1:2 % first two steps of fully implicity
        V = penaltyIteration(V, Vs, Mhat, tol, 0);
    end
    for i = 3:steps % rest with C-N method
        V = penaltyIteration(V, Vs, Mhat, tol, 1/2);
    end
    n = 0;
else
    % variable timestepping
    Tao = dt;
    dTao = dt;
    V_old = V;
    n = 1;
    while(Tao < T)
        if n <=2
            V_new = penaltyIteration(V_old, Vs, M*dTao, tol, 0);
        else
            V_new = penaltyIteration(V_old, Vs, M*dTao, tol, 1/2);
        end
        dTao = varTimeStep(V_new, V_old, dTao, D, dnorm);
        Tao = Tao + dTao;
        V_old = V_new;
        n = n + 1;
    end
    % Now last step to match T
    dTao = T - (Tao  - dTao);
    V_new = penaltyIteration(V_new, Vs, M*dTao, tol, 1/2);
    V = V_new;    
end

% Now interpolate through V to find price based on K
Price = interp1(S,V,K,'nearest');

[delta, gamma] = Greeks(S',V, 2, flag);

end

